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A numerical approach to infrared divergent multi-parton 
phase space integrals 

G. Heinrich'' * 

** II. Institut fiir Thcoretische Physik, Universitat Hamburg, 
Luruper Chaussee 149, 22761 Hamburg, Germany 

It is described how the method of sector decomposition can serve to disentangle overlapping infrared singulari- 
ties, in particular those occurring in the calculation of the real emission part of e^e~ -^ 2 jets and e^e~ —* 3jets 
at NNLO. 



1. INTRODUCTION 

Particle physics nowadays has largely become a 
matter of high precision measurements, which on 
the theory side requires an increasing number of 
loops and legs to be included in the calculations. 

The process e+e" ^jets is particularly inter- 
esting, from an experimental as well as a theoret- 
ical point of view, because of its "clean" initial 
state, such that it can serve for a very accurate 
determination of a^. However, LEP experiments 
already have shown that theoretical predictions at 
next-to-leading order (NLO) are not always suffi- 
cient to match the experimental precision [J , and 
this of course will be even more true for a future 
Linear Collider. 

These facts have triggered a lot of progress 
in the calculation of NNLO corrections in recent 
years ^. The last missing piece for the construc- 
tion of an NNLO Monte Carlo program for the 
process e+e^ -^ 3jets is in fact the double real 
radiation part, which involves phase space inte- 
grations over five partons, where up to two of 
them can become unresolved, leading to infrared 
singularities. 

The conventional way to deal with these singu- 
larities is to establish a subtraction scheme to iso- 
late the divergent part |3I4I5) . The latter then is 
calculated analytically in D = 4 — 2e dimensions, 
leading to 1/e poles which will cancel against the 
ones from the virtual corrections. This proce- 



dure has been applied very successfully in NLO 
calculations. Its generalization to NNLO how- 
ever is far from being straightforward. Neverthe- 
less, subtraction schemes have been proposed in 
the literature 6 7 8 9 lO| , whereas the problem of 
integrating the subtraction terms analytically in 
D dimensions has been solved only for the case 
e+e- ^ 2jets so far |11I1!JI8| . 

What will be suggested here is a new method 
which does not rely on explicit subtraction terms. 
The infrared singularities are isolated in an au- 
tomated way using sector decomposition |13I14| . 
The cancellation of the pole coefficients with the 
ones from the virtual corrections can be verified 
numerically. The method already has been ap- 
plied successfully to the process e^e' -^ 2jets 
|12I15I1(JI17 |. 



2. SECTOR DECOMPOSITION 

The method of sector decomposition acts on 
parameter integrals and serves to factorize singu- 
larities which have an overlapping structure, as 
in the following simple example: 
I — Jq dxi dx2 Xi^~^ [xi + X'2\~^ ■ Decomposing 
the parameter space into two sectors where the 
integration variables are ordered and remapping 
the integration range to the unit square factorizes 
the singularity: 
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[e(a;i-.T2) + e(x2-a;i) 

(1) (2) 



The substitution X2 = Xi t2 in sector (1) and xi = 
X2 ti in sector (2) leads to 

1=1 dxiX^^-'- [ dt2{l+t2)-^ 

Jo Jo 

+ / dx2X-^-' dtit-^-'{l+ti)-K 

Jo Jo 

For more complicated functions, this procedure 
may have to be iterated, but the principle is sim- 
ple and easily automated. This is particularly 
true for multi-loop integrals because they have, 
after Feynman parametrization and integration 
over the loop momenta, the following universal 
form [L is the number of loops, N the number of 
propagators and D the space-time dimension) 

G = {-lfT{N^LDl2) \{dx, (1) 

^ = 1 

^(2.)JV-(L+l)D/2 

4=1 



<5(l-^x. 



T{x,{s,m^})N-LD/2 ' 



where Li and T are polynomials in the Feynman 
parameters and J- also contains kinematic invari- 
ants. Applying the sector decomposition algo- 
rithm ^1] to loop integrals in the form ^ isolates 
the dimensionally regulated poles in terms of fac- 
torizing Feynman parameters. Then subtractions 
of the singularities are carried out, using identi- 
ties like 



/ dxiXi '^J^{xi,x) 
Jo 



dxi T{xi, x) S{xi) + 



/ dxi Xi 
Jo 



[T{xi,x)-T{0,x)] 



(2) 



where x = X2,- ■ ■ ,xn and lima;j^o-?^(a^i7^) is fi- 
nite by construction, such that the second term 
in JSJl is a plus distribution: 
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Doing these subtractions for all Xi results in a 
Laurent series 



I = 



E 

k=-2L 



e'Cfc 



+ 0(6"+^) 



where the order b of expansion in e is in prin- 
ciple only limited by CPU time. However, the 
pole coefficients Ck being sums of complicated 
parameter integrals, their analytical evaluation is 
in general impossible. Therefore they are inte- 
grated numerically. For multi-loop integrals in- 
volving more than one kinematic invariant, Eu- 
clidean points have to be chosen in order to have 
stable numerics. In this way, results have been 
obtained JS| for example for massless 2-loop 4- 
point functions with 2 off-shell legs, where no an- 
alytical results exist yet, all 4-point master inte- 
grals needed for the calculation of 2-loop Bhabha 
scattering with massive fermions (analytical re- 
sults exist for two of them ,19 20 21 ), two-point- 
functions with 4 and 5 loops, and for the pla- 
nar massless 3-loop 4-point-function with on-shell 
legs calculated analytically by V.A. Smirnov {2^ . 

3. PHASE SPACE INTEGRALS 

The phase space integration for the production 
of N massless particles q -^ pi, . . . , p^ can be 
written as 

» Af N 

Y[d^P,S+ip'^)S{q-J2p?) 

j = l 4=1 

^ (2^)A'-I3(JV-1) 21-A' 

•^ j = l ^ 4=1 

At this point one could pick a particular frame 
and integrate over energies Ej and angles 9j , but 
for our purposes it is more convenient to integrate 
over the scaled invariants Sij/q'^, Sij = {pi +pj)'^, 
because in this way the singularities are located 
at the origin of parameter space and no particu- 
lar axis is preferred. The transformation to the 



integration variables 

xi = si2/q^,X2 = 513/9^,2:3 = S2z/q^, 



Xi 



su/q'^,X5 = S24:/q'^,xe = ss^/q'^ 



introduces a Jacobian which is proportional to 
the square root of the determinant of the Gram 
matrix dj — 2piPj. The phase space then takes 
the form^ 



d^l^AT = cW(52)(Ar-l)D/2-^ j J| ^^^. 



i=i 



Us 






e(AAr) (3) 



i=l 

Us = N{N -l)/2 
AN^\detG\iq^r^ 

xV{D- 1)...V{D-N + 1) 
V{D) = 27rT/r(|) . 

3.1. 1-^4 phase space 

As an example, let us consider the integration 
of some squared matrix element IM4P over the 
1 — > 4 partonic phase space, relevant for the cal- 
culation of e+e" -^ 2jets at NNLO: 



d$i^4 

6 



^(4) (^2)313/2-5 



j\{dx,5{l~Y.'^i)\Mi 



'\{xiXq, X2Xz,X'iXi)Y^''^^^ 0(-A) 



(4) 



\{x^y,z)=x +y +z —2{xy + xz + yz) 
The matrix element is of the form 
,2 Vi{x,e) 



iM4r 



(x2 + X4 + xe){x3 + X5 + xe)x4 
, ^2(x,e) ^ 

X2ix2 + X4 + Xq)'^ 



^Note that for N > 6 det G is zero for 4-dimensional mo- 
menta because, after elimination of pg by momentum con- 
servation, the vectors pi to ps will still be linearly depen- 
dent. Therefore we only consider the case N < 6 here. 



where the Vk {x, e) are some polynomials in the 
variables Xi . We again see the sums of Feynman 
parameters in the denominator, corresponding to 
triple invariants Sijk, giving rise to an overlap- 
ping structure. Therefore, the form of the in- 
tegral (0J is very similar to the one in eq. 1^ 
for loop integrals and the overlapping singular- 
ities can be disentangled by the same principle. 
However, there are also very important differ- 
ences to loop integrals. The most important one 
consists in the fact that in phase space integrals, 
non-polynomial structures (square roots) appear. 
For example, solving the constraint — A > 



in Q) for xg leads to 



< xq < xt with 



Xq = {y^X2X^ — ^/xsXi) /xi. The substitution 
xq -^ {xq ~ Xq ) ya + Xq remaps the integration 
range of xg to an integral from to 1 again and 
factorizes the A-term: 



-A] 



-1/2-e 



[xfix^ -xe){xG 



-l/2-£ 



[yeii-ye) 



1-1/2-e 



[a:i(a;([ -Xg)]" 



However, it is possible to eliminate the square 
roots by quadratic transformations, except in fac- 
tors like (1 — yl)"^^^"'^! which do not lead to sin- 
gularities in e and therefore are not subject to 
further sector decomposition. This nice feature 
will be spoiled in the 1^5 case. 

The implementation of sector decomposition 
for the 1^4 phase space served for the calcula- 
tion of all master phase space integrals which are 
needed for any 1 — > 4 process in massless QCD. 
These master integrals have been derived and cal- 
culated analytically as well as numerically in |12| . 

Moreover, the method also can deal with the 
full matrix element without reduction to master 
integrals. This has been demonstrated in [TU). To 
split the calculation into smaller pieces, one can 
write the squared matrix element as a sum over 
different topologies. As the calculation is natu- 
rally parallelized by this subdivision into topolo- 
gies, the overall runtime is given by the most dif- 
ficult topology, which took about 9 hours for a 
precision of 0.1% and less than two hours for a 
precision of 1% on a Pentium IV 2.2 GHz PC. 

3.2. 1^5 phase space 

The 1 — > 5 partonic phase space, relevant for 
the calculation of e+e^ -^ 3jets at NNLO, in- 



volves the integration over 9 independent invari- 
ants: 



d$ 



1^5 
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C, 
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<5(l-^x.)[A5(x)](^-^'/^e(A 



10 
(D-6)/2o, 
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i=l 



c, 



(5) 



V{D - 4) = 27r-7r(-e) 



Note that 

of order e, therefore the integral (0 contains a 
fake singularity in [A5(af)]^ " '' — [A5(x)]~ ~^, 
but this presents no problem for sector decompo- 
sition as the algorithm will extract the singular 
factor and the e-expansion subroutine will take 
the prefactor of order e into account, such that 
the fake singularity will be eliminated automat- 
ically. What is more of a problem are the non- 
polynomial structures which occur here, because 
denominators of the form g{x,y) — a + x + y — 
\/ a? -\- X + y^ where a is a constant, can produce 
a singularity for x,y —^ without having the 
right scaling behaviour amenable to sector de- 
composition. The task is to transform such terms 
away without increasing the complexity of the in- 
tegrand too much. It should be noted that the 
size of the expressions in the 1 — > 5 case is con- 
siderably larger than in the 1 — > 4 case, such that 
it becomes much more important to produce as 
few subsectors as possible. 

The simplest example to calculate is the 5- 
particle phase space volume without any matrix 
element. In |12| a general analytic expression for 
the 1 — > iV phase space volume is given, such 
that the numerical result can be easily checked. 
By sector decomposition, one obtains 

(47r)4'^-7 



d$i- 



r(l-2e)r(2-2e) 



0.00347 



-0.05469e4 
-2.474246^ 



0.44336e^ 
hl0.7283e^ + O(e^ 



which agrees with the analytical result to an accu- 
racy of 0.5% after a runtime of about 10 minutes. 

3.3. One loop plus single real emission 

Apart from the double real emission and the 
two-loop virtual contributions to the cross sec- 
tion of e+e^ -^ jets at NNLO, there is also a 



contribution where one-loop virtual corrections 
are combined with single real emission. In this 
class, the most complicated diagram which can 
occur in the calculation of e"'"e~ — > 2 jets is a box 
graph with one off-shell leg. This type of dia- 
gram can easily be calculated by sector decom- 
position: The one-loop box can be expressed by 
Hypergeometric functions 2Fi{l, — e, 1 — e; Xi/xj). 
Then the parameter representation of the Hyper- 
geometric functions can be used and the result- 
ing one-dimensional parameter integrals can be 
combined with the ones for the 3-particle phase 
space to end up with a 4-dimensional parameter 
integral which can be directly fed into the sector 
decomposition routine. 

For e+e" -^ 3jets, the most complicated one- 
loop diagrams are pentagons with one off-shell 
leg. These could be reduced to boxes by standard 
reduction techniques J23I24| , but as the reduction 
introduces inverse determinants of kinematic ma- 
trices which may lead to numerical instabilities, 
it is more convenient to apply the sector decom- 
position routine for loop integrals directly to the 
pentagon which is of the form 



h 



-r(3 



Y[dz^S{l 



1=1 



5 
i=l 



\r 



■3+e 



-^ = Si2 Z1Z5 + S23 21 (Z3 + Z4 + Z5) 

+ Sl3 Z5{zi + Z2) + Sl4 Zzizi + Z2 + Z3) 
+ S24 Zi{z4 + Z5) + S34 [Zi + Z2){Z4 + Z5) . 

After sector decomposition in the variables Zi , one 
obtains an expression where the poles of the vir- 
tual integral already have been extracted: 

2 

a=0 

pi 4-a 

Po, ^ Y\_'^^iS{U,Sl2,.. .,334) , 

•^0 t=l 

hm g ^0 . 

This expression can then be inserted into the 4- 
particle phase space and one can proceed with 
decomposition in the scaled invariants Xi, . . . ,Xq. 
Note that no problems with thresholds will occur 
here as the kinematics is such that all invariants 
Sij are non-negative. 



4. SUMMARY AND OUTLOOK 

The automated sector decomposition algorithm 
is a powerful method to isolate overlapping in- 
frared poles and to calculate numerically not only 
multi-loop integrals, but also phase space inte- 
grals where some of the particles can become 
theoretically unresolved, leading to infrared sin- 
gularities. In particular, the method allows the 
calculation of the one-loop plus single real emis- 
sion and the double real emission contribution 
to e^e~ ^ 2 or 3 jets at NNLO without hav- 
ing to establish a subtraction scheme and to in- 
tegrate analytically over complicated subtraction 
terms. The inclusion of a measurement func- 
tion also does not present a problem, as has been 
demonstrated already in \n\ , such that a fully dif- 
ferential Monte Carlo program can be constructed 
based on this method. The only drawback of the 
method is the fact that it generates a large num- 
ber of functions, but it has been shown already 
that in the case of e+e" -^ 2 jets at NNLO, 
this does not lead to unacceptable integration 
times. Further, the functions are numerically 
well-behaved by construction. How the NNLO 
calculation of the process e^e~ — > 3 jets with 
this method performs numerically will turn out 
in the near future. 

The generalization to other processes than 
e^e~ annihilation is feasible, but cases where 
some of the kinematic invariants take negative 
values cannot be treated without further deve- 
lopment of the method. 
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